###########################################################
#       BIAS WILL FIND A WAY	  	 		  #
#	replication script				  #
#	Author: Martin Bisgaard 			  #
#	Date: 02-03-2015
###########################################################

rm(list=ls())

##PREDICT FUNCTION WITH ALTERNATE SEs ##
# function is used to check robustness for specificatoins w/o robust SEs
predict.new<- function(x,clcov,newdata){
	if(missing(newdata)){ newdata <- x$model }
	tt <- terms(x)
	Terms <- delete.response(tt)
	m.mat <- model.matrix(Terms,data=newdata)
	m.coef <- x$coef
	fit <- as.vector(m.mat %*% x$coef)
	se.fit <- sqrt(diag(m.mat%*%clcov%*%t(m.mat)))
	return(list(fit=fit,se.fit=se.fit))
	}



require(TTR)
require(sandwich)
require(memisc)
require(foreign)
require(car)
require(stargazer)
require(lmtest)


load("newdat.RData")
dim(newdat)#should be 87360 by 30

###############################
#	CODING VARIABLES		#


newdat$readnwsp<-recode(newdat$readnwsp,'8=NA')
newdat$education<-recode(newdat$education,'9=NA')
newdat$income<-recode(newdat$income,'11:15=NA;98=NA')

newdat$rpecon2<-recode(newdat$rpecon,'1=1;2=0;3=NA;4=0;5=0;8=NA;9=NA')
newdat$rpecon3<-recode(newdat$rpecon,'1=1;2=NA;3=NA;4=0;5=0;8=NA;9=NA')
newdat$rpecon4<-recode(newdat$rpecon,'1=1;2=0;3=1;4=0;5=0;8=NA;9=NA')
newdat$rpecon<-recode(newdat$rpecon,'1=1;2=0;3=0;4=0;5=0;8=NA;9=NA')

newdat$econprob<-recode(newdat$econprob,'12=NA');newdat$econprob<-newdat$econprob/10
newdat$rpgov<-recode(newdat$rpgov,'9=NA')
newdat$coneu<-recode(newdat$coneu,'8=NA;9=NA')

newdat$unemp<-recode(newdat$unemp,'3:4=1;else=0;99=NA')
newdat$gender<-recode(newdat$gender,'1=0;2=1;8=NA')
newdat$mortgage<-recode(newdat$mortgage,'1=0;2=1')

newdat$party.str.full<-newdat$party.str
newdat$party.str<-recode(newdat$party.str,'1=3;3=1;4=NA;8=NA');newdat$party.str<-(newdat$party.str-1)/2
newdat$coneu<-(newdat$coneu-1)/3

newdat$polatt<-recode(newdat$polatt,'12=NA');newdat$polatt<-newdat$polatt/10


newdat$natretro<-recode(newdat$natretro,'1=5;2=4;3=3;4=2;5=1;6=3;8=NA;9=NA')
newdat$natretro <- (newdat$natretro-1)/4;table(newdat$natretro)

newdat$consgov<-0
newdat$consgov[newdat$date %in% 72:97]<-1
table(newdat$consgov)


###################
#	PARTY ID	#

newdat$party.id2<-NA
newdat$party.id2[newdat$party.id==2]<-0
newdat$party.id2[newdat$party.id==1]<-1

newdat$party.id4<-NA
newdat$party.id4[newdat$party.id==2]<-0 #lab
newdat$party.id4[newdat$party.id==1]<-1 #cons
newdat$party.id4[newdat$party.id==3]<-2 #libdem
newdat$party.id4[newdat$party.id %in% c(10)]<-3 #no-one + D/K=11
newdat$party.id4<-as.factor(newdat$party.id4)
table(newdat$party.id4)

newdat$party.id5<-NA
newdat$party.id5[newdat$party.id==2]<-0 #lab
newdat$party.id5[newdat$party.id==1]<-1 #cons
newdat$party.id5[newdat$party.id==3]<-2 #libdem
newdat$party.id5[newdat$party.id %in% c(4,6,7,8,9)]<-3 #Other
newdat$party.id5[newdat$party.id %in% c(10,11)]<-4 #None + DK
newdat$party.id5<-as.factor(newdat$party.id5)
table(newdat$party.id5)

newdat$vote2<-recode(newdat$vote2,'2=0;1=1;else=NA');table(newdat$vote2)
newdat$primesat<-recode(newdat$primesat,'1=4;2=3;3=2;4=1;8=NA');table(newdat$primesat)
newdat$primesat<-(newdat$primesat-1)/3;table(newdat$primesat)


newdat$age<-recode(newdat$age,'98:99=NA')
newdat$age<-newdat$age-1 #normalize so 0 is 1900




#########################
# 	ANALYSIS	#
#########################


#	FIGURE 2:				#
# 	TIME-VARYING DYNAMICS OF NATRETRO 	#


lmfit<-lm(natretro~party.id2*date+party.str+polatt+gender+age+income+mortgage+unemp+education+readnwsp,data=newdat)
summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T),
	"readnwsp"=median(newdat$readnwsp,na.rm=T)
))



#vmat<-vcovHC(lmfit,type="HC1") # ROBUST SEs do not make any substantive difference for the results (uncomment and set vmat as the Robust var-covar to check)

vmat<-vcov(lmfit) #conventional var-covar estimator is used


preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))

coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

################
# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")


############
# PLOTTING # 


coefs.out<-rbind(coefs.out[1:13,],NA,coefs.out[14:28,],NA,coefs.out[29:71,]) #adding blanks for the two missing months


col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)



pdf('retro.pdf',width=5.75,height=6.25)

par(mfrow=c(1,1))
par(omd=c(0,.95,.375,1))

plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Predicted Mean",yaxt="n",
	col=col1)
axis(2,at=seq(0,1,.25),las=2)
axis(4,at=seq(0,1,.25),labels=c("Got a lot","   Got"," Around","   Got","Got a lot"),cex.axis=.75,las=2,padj=-.25,hadj=.1)
axis(4,at=seq(0,1,.25),labels=c("  worse","  worse","the same","  better","  better"),cex.axis=.75,las=2,padj=.75,hadj=.1,tick=F)

segments(1:73,coefs.out[,1],1:73,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:73,coefs.out[,1],1:73,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:73,coefs.out[,3],1:73,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:73,coefs.out[,3],1:73,coefs.out[,3]-1.97*coefs.out[,4],col=col1)

points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)


lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,1]),n=3))
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,3]),n=3))


abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=2)

legend(-3,1.05,pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,.95,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,.95,.375,1),new=TRUE)

text(69,.9,labels="Northern Rock",font=1,cex=1,pos=2)
arrows(46,.88,42,.8,length=.075)
#legend("bottomleft",lty=c(1,1),pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour identifier","Conservative Identifier"),cex=.9,bty="n")


par(omd=c(0,.95,0,.625),new=TRUE)

plot(coefs.out[,5],ylim=c(-.1,.55),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Marginal Effect of Party Id.",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),las=2)
abline(h=0,lty=2)


abline(h=.5)
mtext(3,text="Partisan Differences",line=-1,font=2)


par(omd=c(0,.95,0,.585),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,.95,0,.625),new=TRUE)



segments(1:73,coefs.out[,5],1:73,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:73,coefs.out[,5],1:73,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,5]),n=2))
points(coefs.out[,5],pch=21,bg=col2,col=col2)

axis(1,c(10,22,34,46,58,70),labels=c(2005,"",2007,"",2009,""),cex.axis=1)


dev.off()



#########################################################
#		FIGURE 3				#
# 	TIME-VARYING DYNAMICS OF ECONOMIC PERCEPTIONS 	#
#SUBSET APPROACH TO PARTY.STR interaction		#
#							#


#WEAK IDENTIFIERS##

lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education+readnwsp
	,data=newdat[newdat$party.str==0,])
summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T),
	"readnwsp"=median(newdat$readnwsp,na.rm=T)
))


vmat<-vcov(lmfit)



preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

################
# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")

coefs.out<-rbind(coefs.out[1:13,],NA,coefs.out[14:28,],NA,coefs.out[29:71,]) #adding blanks for the two missing months (for plotting)
coefsWEAK<-coefs.out 	#storing





######MIDDLING IDENTIFIERS########



lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education+readnwsp
	,data=newdat[newdat$party.str==0.5,])
summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T),
	"readnwsp"=median(newdat$readnwsp,na.rm=T)
))

vmat<-vcov(lmfit)

preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

################
# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)
coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")


coefs.out<-rbind(coefs.out[1:13,],NA,coefs.out[14:28,],NA,coefs.out[29:71,]) #adding blanks for the two missing months (for plotting)
coefsMID<-coefs.out





# 	STRONG IDENTIFIERS #


lmfit<-lm(natretro~party.id2*date+polatt+gender+age+income+mortgage+unemp+education+readnwsp
	,data=newdat[newdat$party.str==1,])
#summary(lmfit)




#####################
# CONDITIONAL MEANS #

preddat<-data.frame(expand.grid("party.id2"=c(0,1),"date"=as.factor(1:71),
	"party.str"=median(newdat$party.str,na.rm=T),
	"polatt"=median(newdat$polatt,na.rm=T),
	"gender"=median(newdat$gender,na.rm=T),
	"age"=mean(newdat$age,na.rm=T),
	"income"=median(newdat$income,na.rm=T),
	"mortgage"=median(newdat$mortgage,na.rm=T),
	"unemp"=median(newdat$unemp,na.rm=T),
	"education"=median(newdat$education,na.rm=T),
	"readnwsp"=median(newdat$readnwsp,na.rm=T)
))

vmat<-vcov(lmfit)

preddat.out<-cbind(preddat,predict.new(lmfit,vmat,preddat))


coefs.out<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting
coefs.out[,1]<-preddat.out$fit[preddat.out$party.id2==0]
coefs.out[,2]<-preddat.out$se.fit[preddat.out$party.id2==0]
coefs.out[,3]<-preddat.out$fit[preddat.out$party.id2==1]
coefs.out[,4]<-preddat.out$se.fit[preddat.out$party.id2==1]

################
# PARTISAN GAP #

est<-matrix(c(coef(lmfit)["party.id2"],
	rep(coef(lmfit)["party.id2"],70)+
	coef(lmfit)[grep("^party.id2:date",names(coef(lmfit)))]
		))


###########################
# SEs for marginal effect #


b1<-rep(vmat["party.id2","party.id2"],70)
b3<-diag(vmat)[grep("^party.id2:date",rownames(vmat))]
b1b3<-vmat["party.id2",grep("^party.id2:date",rownames(vmat))] #extract covar

se<-cbind(b1,b3,b1b3);rownames(se)<-rep("",70) #arrange in matrix form

se.out<-apply(se,1,function(x) sqrt(x[1]+x[2]+2*x[3]))
se.out<-c(b1[1],se.out)

################
# Final output #

est.out<-cbind(est,se.out)
coefs.out[,5:6]<-est.out;coefs.out<-round(coefs.out,digits=5)

coefs.out[,5]<- -coefs.out[,5]
colnames(coefs.out)<-c("LAB","SE","CON","SE","GAP","SE")

coefs.out<-rbind(coefs.out[1:13,],NA,coefs.out[14:28,],NA,coefs.out[29:71,]) #adding blanks for the two missing months (for plotting)
coefsSTRONG<-coefs.out




#############
# PLOTTING	#

col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)



pdf('natretro-strengthcomposite.pdf',width=10,height=6)


d1<-.4-0.05/2
d2<-.3-0.05/2
d3<-.6-0.05/2
d4<-.7-0.05/2

td<-1.75
subd<-0.25


par(mfrow=c(1,1))
par(omd=c(0,d1,.375,1))

coefs.out<-coefsWEAK

plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Predicted Mean",yaxt="n",xaxt="n",
	col=col1)
axis(1,c(10,22,34,46,58,70),labels=NA)
axis(2,at=seq(0,1,.25),las=2)

mtext(3,text='A',line=td,font=2,cex=1.5)
mtext(3,text='Weak Identifiers (N=8,933) ',line=subd,font=1,cex=1)


segments(1:73,coefs.out[,1],1:73,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:73,coefs.out[,1],1:73,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:73,coefs.out[,3],1:73,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:73,coefs.out[,3],1:73,coefs.out[,3]-1.97*coefs.out[,4],col=col1)

points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)

lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,1]),n=3))
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,3]),n=3))


abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)

legend(-3,1.05,pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour Identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,d1,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,d1,.375,1),new=TRUE)



par(omd=c(0,d1,0,.625),new=TRUE)

plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Marginal Effect of Party Id.",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),las=2)
abline(h=0,lty=2)


abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)


par(omd=c(0,d1,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,d1,0,.625),new=TRUE)



segments(1:73,coefs.out[,5],1:73,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:73,coefs.out[,5],1:73,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,5]),n=3))
points(coefs.out[,5],pch=21,bg=col2,col=col2)

axis(1,c(10,22,34,46,58,70),labels=c(2005,"",2007,"",2009,""),cex.axis=1)


##ADDING MIDDLING CATEGORY###

par(omd=c(d2,d4,.375,1),new=T)

coefs.out<-coefsMID

plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1)
axis(1,c(10,22,34,46,58,70),labels=NA)
axis(2,labels=NA)

mtext(3,text='B',line=td,font=2,cex=1.5)
mtext(3,text='Middling Identifiers (N=17,437)',line=subd,font=1,cex=1)


segments(1:73,coefs.out[,1],1:73,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:73,coefs.out[,1],1:73,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:73,coefs.out[,3],1:73,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:73,coefs.out[,3],1:73,coefs.out[,3]-1.97*coefs.out[,4],col=col1)

points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)

lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,1]),n=3))
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,3]),n=3))


abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)

text(0,.9,labels="Northern Rock",font=1,cex=.8,pos=4)
arrows(29,.875,42,.8,length=.075)



par(omd=c(d2,d4,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d2,d4,.375,1),new=TRUE)



par(omd=c(d2,d4,0,.625),new=TRUE)

plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)
abline(h=0,lty=2)


abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)


par(omd=c(d2,d4,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d2,d4,0,.625),new=TRUE)



segments(1:73,coefs.out[,5],1:73,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:73,coefs.out[,5],1:73,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,5]),n=3))
points(coefs.out[,5],pch=21,bg=col2,col=col2)

axis(1,c(10,22,34,46,58,70),labels=c(2005,"",2007,"",2009,""),cex.axis=1)




##ADDING STRONG CATEGORY###


par(omd=c(d3,.95,.375,1),new=T)

coefs.out<-coefsSTRONG

plot(coefs.out[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1)
axis(1,c(10,22,34,46,58,70),labels=NA)
axis(2,labels=NA)

axis(4,at=seq(0,1,.25),labels=c("Got a lot","   Got"," Around","   Got","Got a lot"),cex.axis=.75,las=2,padj=-.25,hadj=.1,xpd=T)
axis(4,at=seq(0,1,.25),labels=c("  worse","  worse","the same","  better","  better"),cex.axis=.75,las=2,padj=.75,hadj=.1,tick=F,xpd=T)


mtext(3,text='C',line=td,font=2,cex=1.5)
mtext(3,text='Strong Identifiers (N=6,522)',line=subd,font=1,cex=1)


segments(1:73,coefs.out[,1],1:73,coefs.out[,1]+1.97*coefs.out[,2],col=col1);segments(1:73,coefs.out[,1],1:73,coefs.out[,1]-1.97*coefs.out[,2],col=col1)
segments(1:73,coefs.out[,3],1:73,coefs.out[,3]+1.97*coefs.out[,4],col=col1);segments(1:73,coefs.out[,3],1:73,coefs.out[,3]-1.97*coefs.out[,4],col=col1)

points(coefs.out[,1],pch=21,bg=col1,col=col1)
points(coefs.out[,3],pch=21,col=col1)

lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,1]),n=3))
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,3]),n=3))


abline(h=1.05)
mtext(3,text="Retrospective Economic Perceptions",line=-1,font=1,cex=.8)


par(omd=c(d3,.95,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d3,.95,0,.625),new=TRUE)

plot(coefs.out[,5],ylim=c(-.1,.6),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt='n',
	col="white")
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)
abline(h=0,lty=2)


abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)


par(omd=c(d3,.95,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d3,.95,0,.625),new=TRUE)



segments(1:73,coefs.out[,5],1:73,coefs.out[,5]+1.97*coefs.out[,6],col=col2);segments(1:73,coefs.out[,5],1:73,coefs.out[,5]-1.97*coefs.out[,6],col=col2)
lines(c(1:13,15:29,31:73),SMA(na.omit(coefs.out[,5]),n=3))
points(coefs.out[,5],pch=21,bg=col2,col=col2)

axis(1,c(10,22,34,46,58,70),labels=c(2005,"",2007,"",2009,""),cex.axis=1)

dev.off()





#############################################################
# TIME-VARYING DYNAMICS OF ATTRIBUTIONS OF RESPONSIBILITY   #
#		FIGURE 4				    #


coefs.out<-matrix(NA,nc=6,nr=97)
boot<-1000 #set number of reps for AME, WARNING: runs slow with 1000
set.seed(123)

lmfit<-glm(rpecon2~party.id2*date+party.str+polatt+coneu+gender+age+income+mortgage+unemp+education+readnwsp,data=newdat,
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=71) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","readnwsp","mortgage","education","income","gender","age","polatt","unemp","party.id2","date")

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #extract var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty

	#set up date frame with observed values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:71){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$date==i,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$date==i,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp




###################
#	Plotting	#

coefsrp[,5]<- -coefsrp[,5] #reverse gap coefficient

coefsrp<-rbind(coefsrp[1:13,],NA,coefsrp[14:28,],NA,coefsrp[29:71,]) #adding blanks for the two missing months


col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)

pdf('rpecon.pdf',width=5.75,height=6.25)

par(mfrow=c(1,1))
par(omd=c(0,1,.375,1))

plot(coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Predicted Probability",
	col=col1)

segments(1:73,coefsrp[,1],1:73,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(1:73,coefsrp[,1],1:73,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(1:73,coefsrp[,3],1:73,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(1:73,coefsrp[,3],1:73,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(coefsrp[,1],pch=21,bg=col1,col=col1)
points(coefsrp[,3],pch=21,col=col1)

lines(c(1:13,15:29,31:73),SMA(na.omit(coefsrp[,1]),n=3))
lines(c(1:13,15:29,31:73),SMA(na.omit(coefsrp[,3]),n=3))


abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=2)

legend("bottomleft",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,1,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,1,.375,1),new=TRUE)

text(35,.975,labels="Northern Rock",font=1,cex=1,pos=2)
arrows(34,.975,42,.95,length=.075)

par(omd=c(0,1,0,.625),new=TRUE)

plot(coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Marginal Effect",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=2)

segments(1:73,coefsrp[,5],1:73,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(1:73,coefsrp[,5],1:73,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(c(1:13,15:29,31:73),SMA(na.omit(coefsrp[,5]),n=3))
points(coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(0,1,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,1,0,.625),new=TRUE)

axis(1,c(10,22,34,46,58,70),labels=c(2005,"",2007,"",2009,""),cex.axis=1)



dev.off()
 




###############################################
#	FIGURE 5			      #
# Attribution conditional on strength of PID  #
# BASED ON 6month date variable		      #

set.seed(123)
boot<-1000 #set number of reps for AME and APP

#collapse date into 6month bins
#note that two of the bins include 5 waves (and not 6) because the CMS was not conducted here. 
newdat$datehy<-recode(as.numeric(newdat$date),
	'1:6=0;
	7:12=1;
	13:17=2; 
	18:23=3;
	24:28=4; 
	29:34=5;
	35:40=6;
	41:46=7;
	47:52=8;
	53:58=9;
	59:64=10;
	65:71=11;else=NA')
table(newdat$date,newdat$datehy)#check joint distribution with date
newdat$datehy<-factor(newdat$datehy)


##############
#  WEAK ID   #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education+readnwsp,
	data=newdat[newdat$party.str==0,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting

varnam<-c("coneu","party.str","readnwsp","mortgage","education","income","gender","age","polatt","unemp","party.id2","datehy")

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x)
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpWEAK<-coefsrp

#write.csv(coefsrp,"coefsrpecon2FULL.csv")
#write.csv(coefsrp[65:11,],"coefsrpecon2GOV.csv")





##############
#  MIDDLING ID   #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education+readnwsp,
	data=newdat[newdat$party.str==.5,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting


################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #cluster robust var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpMID<-coefsrp

#write.csv(coefsrp,"coefsrpecon2FULL.csv")
#write.csv(coefsrp[65:11,],"coefsrpecon2GOV.csv")






################
#  STRONG ID   #

lmfit<-glm(rpecon2~party.id2*datehy+polatt+coneu+gender+age+income+mortgage+unemp+education+readnwsp,
	data=newdat[newdat$party.str==1,],
	family=binomial(link="probit"))

mtable(lmfit)

coefsrp<-matrix(NA,nc=6,nr=12) #matrix to hold overall output for plotting

	vmat<-cluster.vcov(x,dat$date[rownames(newdat)%in%rownames(x$model)]) #cluster robust var-covar matrix
	coeftest(x,vcov=vmat)

################
# APP and AME  #

	x<-lmfit 	#Set model
	vmat<-vcov(x) #cluster robust var-covar matrix
	samp<-mvrnorm(n=boot,mu=coef(x),Sigma=vmat) #simulate estimation uncertainty
dim(vmat)
	#set up datehy frame with observered values
	datx1<-newdat[complete.cases(newdat[,varnam]),varnam];datx1$party.id2<-1 #data where all x's are set to 1
	datx0<-newdat[complete.cases(newdat[,varnam]),varnam];datx0$party.id2<-0 #set to 0

	ame.samp<-rep(NA,boot)
	app1<-rep(NA,boot)
	app0<-rep(NA,boot)
	ame<-rep(NA,boot)

	for (i in 1:12){
	for (k in 1:boot){
		x$coefficients<-samp[k,]
		app0[k]<-mean(predict(x,newdat=datx0[datx0$datehy==i-1,],type="response"))
		app1[k]<-mean(predict(x,newdat=datx1[datx1$datehy==i-1,],type="response"))
		ame[k]<-app0[k]-app1[k]
	}
		coefsrp[i,]<-c(mean(app0),sd(app0),mean(app1),sd(app1),mean(ame),sd(ame)) #storing time varying APPs and AMEs
	}

coefsrp
coefsrpSTRONG<-coefsrp

#write.csv(coefsrp,"coefsrpecon2FULL.csv")
#write.csv(coefsrp[65:11,],"coefsrpecon2GOV.csv")



#########################
#	Plotting 		#

coefsrpWEAK[,5]<- -coefsrpWEAK[,5] #reverse partisan gap coefficient
coefsrpMID[,5]<- -coefsrpMID[,5]
coefsrpSTRONG[,5]<- -coefsrpSTRONG[,5] 




col1<-rgb(.2,.2,.2,alpha=.8)
col2<-rgb(.4,.4,.4,alpha=.8)
pdf('rpecon-strengthcomposite.pdf',width=10,height=6.5)


d1<-.4
d2<-.3
d3<-.6
d4<-.7

td<-1.75
subd<-0.25

par(mfrow=c(1,1))
par(omd=c(0,d1,.375,1))

coefsrp<-coefsrpWEAK


xval<-c(seq(6,71,6)-2.5,69.5)
adjtick<-c(10,22,34,46,58,70,82,94)

plot(xval,coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Predicted Probability",
	col=col1)

mtext(3,text='A',line=td,font=2,cex=1.5)
mtext(3,text='Weak Identifiers (N=4,708)',line=subd,font=1,cex=1)

segments(xval,coefsrp[,1],xval,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(xval,coefsrp[,1],xval,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(xval,coefsrp[,3],xval,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(xval,coefsrp[,3],xval,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(xval,coefsrp[,1],pch=21,bg=col1,col=col1)
points(xval,coefsrp[,3],pch=21,col=col1)

lines(xval,coefsrp[,1])
lines(xval,coefsrp[,3],lty=1)

axis(1,adjtick,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)

legend("bottomleft",pch=c(21,1),pt.bg=c(col1,NA),col=col1,legend=c("Labour Identifier","Conservative Identifier"),cex=.9,bty="n")

par(omd=c(0,d1,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,d1,.375,1),new=TRUE)

text(36,.975,labels="Northern Rock",font=1,cex=.8,pos=2,xpd=T)
arrows(35,.925,42,.8,length=.075)



par(omd=c(0,d1,0,.625),new=TRUE)

plot(xval,coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="Avg. Marginal Effect",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(xval,coefsrp[,5],xval,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(xval,coefsrp[,5],xval,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(xval,coefsrp[,5])
points(xval,coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(0,d1,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(0,d1,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)


########MIDDLING IDers#######
coefsrp<-coefsrpMID

par(omd=c(d2,d4,.375,1),new=TRUE)
plot(xval,coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1,add=TRUE)

mtext(3,text='B',line=td,font=2,cex=1.5)
mtext(3,text='Middling Identifiers (N=9,875)',line=subd,font=1,cex=1)

segments(xval,coefsrp[,1],xval,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(xval,coefsrp[,1],xval,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(xval,coefsrp[,3],xval,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(xval,coefsrp[,3],xval,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(xval,coefsrp[,1],pch=21,bg=col1,col=col1)
points(xval,coefsrp[,3],pch=21,col=col1)

lines(xval,coefsrp[,1])
lines(xval,coefsrp[,3],lty=1)

axis(1,adjtick,labels=NA)
axis(2,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)


par(omd=c(d2,d4,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d2,d4,.375,1),new=TRUE)


par(omd=c(d2,d4,0,.625),new=TRUE)

plot(xval,coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(xval,coefsrp[,5],xval,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(xval,coefsrp[,5],xval,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(xval,coefsrp[,5])
points(xval,coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(d2,d4,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d2,d4,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)



####### 	ADDING STRONG ID GRAPH #############

coefsrp<-coefsrpSTRONG

par(omd=c(d3,1,.375,1),new=TRUE)
plot(xval,coefsrp[,1],ylim=c(0,1.15),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col1,add=TRUE)

mtext(3,text='C',line=td,font=2,cex=1.5)
mtext(3,text='Strong Identifiers (N=3,920)',line=subd,font=1,cex=1)

segments(xval,coefsrp[,1],xval,coefsrp[,1]+1.97*coefsrp[,2],col=col1);segments(xval,coefsrp[,1],xval,coefsrp[,1]-1.97*coefsrp[,2],col=col1)
segments(xval,coefsrp[,3],xval,coefsrp[,3]+1.97*coefsrp[,4],col=col1);segments(xval,coefsrp[,3],xval,coefsrp[,3]-1.97*coefsrp[,4],col=col1)

points(xval,coefsrp[,1],pch=21,bg=col1,col=col1)
points(xval,coefsrp[,3],pch=21,col=col1)

lines(xval,coefsrp[,1])
lines(xval,coefsrp[,3],lty=1)

axis(1,adjtick,labels=NA)
axis(2,labels=NA)

abline(h=1.05)
mtext(3,text="Attribution of Responsibility",line=-1,font=1,cex=.8)


par(omd=c(d3,1,.375,.95),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d3,1,.375,1),new=TRUE)


par(omd=c(d3,1,0,.625),new=TRUE)

plot(xval,coefsrp[,5],ylim=c(-.5,.65),xaxt="n",type="p",lwd=1,xlab="",las=2,
	ylab="",yaxt="n",
	col=col2)
axis(2,at=seq(-.4,.4,.2),labels=NA,las=2)

abline(h=.55)
mtext(3,text="Partisan Differences",line=-1,font=1,cex=.8)

segments(xval,coefsrp[,5],xval,coefsrp[,5]+1.97*coefsrp[,6],col=col2);segments(xval,coefsrp[,5],xval,coefsrp[,5]-1.97*coefsrp[,6],col=col2)
lines(xval,coefsrp[,5])
points(xval,coefsrp[,5],pch=21,bg=col2,col=col2)
abline(h=0,lty=2)

par(omd=c(d3,1,0,.575),new=TRUE)
abline(v=42,lty=2)
par(omd=c(d3,1,0,.625),new=TRUE)


axis(1,adjtick,labels=c(2005,"",2007,"",2009,"",2011,""),cex.axis=1)

dev.off()




